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It was recently shown that the exact factorization of the electron-nuclear wavefunction allows the 
construction of a Schrodinger equation for the electronic system, in which the potential contains ex¬ 
actly the effect of coupling to the nuclear degrees of freedom and any external fields. Here we study 
the exact potential acting on the electron in charge-resonance enhanced ionization in a model one¬ 
dimensional molecule. We show there can be significant differences between the exact potential 
and that used in the traditional quasistatic analyses, arising from non-adiabatic coupling to the nu¬ 
clear system, and that these are crucial to include for accurate simulations of time-resolved ionization 
dynamics and predictions of the ionization yield. 


Ionization is a fundamental process in the strong-field 
physics of atoms and molecules, lying at the heart of 
many fascinating phenomena such as high harmonic 
generation. Coulomb explosion, laser induced electron 
diffraction, and molecular orbital tomography. The ion¬ 
ization rate from a molecule can be several orders of 
magnitude higher than the rate from the constituent 
atoms at a critical range of internuclear separations. This 
phenomenon, termed charge-resonance enhanced ion¬ 
ization (CREI), was first theoretically predicted about 
twenty years ago QHU, and soon after verified experi¬ 
mentally 15HZI- 

The enhancement in the ionization rate has been 
explained by a quasistatic argument in the pioneer¬ 
ing works of Ref (2} [4} IStiTTl, treating the nuclei as 
instantaneously-fixed point particles, with the electrons 
following the combined potential from the laser field 
and the electrostatic attraction of the nuclei. In any ex¬ 
periment however, the nuclei are neither frozen, nor are 
they point particles; instead their motion can be strongly 
coupled to the electron dynamics and accounting for the 
coupled electron-ion quantum dynamics can be essen¬ 
tial m. Further, the electron does not simply follow 
the field adiabatically, evidenced by the multiple subcy¬ 
cle ionization bursts recently revealed in Refs. 1 131 ITU . 
A handful of calculations treating the full dynamics of 
quantum nuclei and the electron in H 2 f j4| [8|, as well 
as a few experiments, have verified that the essential 
CREI phenomenon remains robust, although the de¬ 
tails are altered by the ionic dynamics. For example, in 
Refs. 13 it was found that nuclear motion washes out 
the two-peak structure predicted in the frozen-nuclei 
analysis of CREI 13 into a single broad peak. Further, 
the CREI enhancement is subdued if, during the ex¬ 
periment, only little of the nuclear density reaches the 
critical internuclear separation |T5j], Hence, to properly 
understand, model, and predict the experiment, a fully 
time-dependent picture of coupled electronic and ionic 
motion is needed. 


In this Letter, we utilize the exact factorization ap¬ 
proach lUSUm to investigate the dynamics of the elec¬ 
tron during the CREI process, as influenced by the field 
and nuclear wavepacket dynamics. In particular, we 
study the exact time-dependent potential that drives the 
electron and fully accounts for coupling to both the field 
and the dynamical nuclei. This concept has been intro¬ 
duced by the exact factorization in its inverse form 111 81 , 
and here we show that, for the CREI process, this ex¬ 
act potential can be remarkably different from the qua¬ 
sistatic potential, or even from a modified quasistatic 
potential that accounts for the width and splitting of the 
nuclear wavepacket. This indicates the need for dynam¬ 
ical electron-nuclear correlation effects to be included 
in the calculation. Further, we identify a new measure 
of ionization appropriate for fully dynamical studies, 
which indicates the regions of the nuclear wavepacket 
associated with the ionizing electron. 

By restricting the motion of the nuclei and the elec¬ 
tron in the Elf molecule to the polarization direction of 
the laser field the problem can be modeled with a one¬ 
dimensional Hamiltonian featuring "soft-Coulomb" in¬ 
teractions |20| (atomic units, e = m = h = 1, are used 
throughout the article, unless otherwise noted): 
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where R and z are the internuclear distance and the elec¬ 
tronic coordinate as measured from the nuclear center- 
of-mass, respectively. The proton mass is denoted as M 
while y e = (2M)/(2M + 1) is the electronic reduced 
mass. The laser field, within the dipole approximation, 
is represented by Vi(z, t) = q e zE(t) where E{t) denotes 
the electric field amplitude and the reduced charge q e = 
(2 M + 2)/(2 M + 1). Such reduced-dimensional models 
have been shown to qualitatively reproduce experimen- 








FIG. 1. Upper panel: Ionization probability, dissociation prob¬ 
ability , average internuclear distance ,{R), as a function of 
number of cycles t/T. The inset depicts the ionization rate 
for: clamped-nucleus (CN) calculation, as a function of the in¬ 
ternuclear distance (black dashed curve), the exact molecular 
TDSE as a function of (R) (red full curve), and the H atom 
(dash dotted line). Middle panel: contour plot of the time- 
dependent nuclear density. Lower panel: contour plot of time- 
dependent electronic density. 


tal results (see Ref. m for example). Here we first study 
the dynamics of the system subject to a 50-cycle pulse of 
wavelength A = 800 nm (u> = 0.0569 a.u.) and intensity 
I = 2 x 10 14 W/cm 2 , with a sine-squared pulse envelope. 

Setting the ground-state of the molecule as the ini¬ 
tial state, we first solve the time-dependent Schrodinger 
equation (TDSE) for Hamiltonian in Eq. (Jljl numerically 
exactly. The upper panel of Figure [lj shows the ioniza¬ 
tion probability, dissociation probability and the aver¬ 
age internuclear distance (R) as a function of number 
of optical cycles (t/T) 1221 . Here T denotes duration 
of one cycle which in this case is T = 2.67 fs. We see 
from the figure that ionization is rapidly onset as we 
approach the middle of the pulse, slowing down later 
while the field decreases. The nuclei dissociate primar¬ 
ily via Coulomb explosion following ionization. Note 
that the majority of the ionization occurs when (R) is 
between 4-5.5 a.u. The inset of the upper panel of Fig¬ 
ure (JTJ shows the ionization rate which in the fully dy¬ 
namical calculation, is plotted against the average inter¬ 
nuclear separation ( R)(t ), while in the clamped-nuclei 
(CN) case, it is calculated and plotted for each fixed in¬ 
ternuclear separation R 11231 . In the CN case, the peak 
near 6.5 a.u. is usually identified with CREI while that 
near 5 a.u. is associated with a symmetry-breaking elec¬ 
tron localization [ 2 ) • The exact ionization rate, however, 
has a single broad peak centered between 4 a.u. and 5 
a.u., and is smaller than that of the CN calculations, but 
still higher than the atomic rate, similar to the observa¬ 
tions in Refs. niia. 

One should keep in mind that representing the ioniza¬ 
tion rate in the fully-dynamical case as a function of just 


the average internuclear separation has however only 
limited meaning: the nuclear charge distribution, plot¬ 
ted in the middle panel of Figure]!] bifurcates, so consid¬ 
ering ionization simply as a function of the average sep¬ 
aration does not properly indicate the internuclear sepa¬ 
rations at which the ionization rate is enhanced. A large 
fragment of the nuclear density remains localized, oscil¬ 
lating around the equilibrium separation, while another 
part begins to dissociate, soon after the ionization is on¬ 
set. This can be seen clearly in the plot of the electronic 
density plotted in the lower panel of Figure [l] Hence a 
dynamical picture of the CREI that accounts for the cou¬ 
pling to the nuclear distribution as it changes in time is 
desirable. 

A complete picture of the electronic dynamics cou¬ 
pled to the non-classical nuclei is provided within 
the exact factorization framework in its inverse for¬ 
mulation 118 |: Complementary to the direct factoriza¬ 
tion of Refs. IBS CEE the exact electron-nuclear wave- 
function 41 (r, R, t) that solves the full electron-nuclear 
TDSE can be exactly written as a product 't (r, R, t) = 
$(r, f)xr(R,t)/ where ‘l>(r, t) may be interpreted as the 
electronic wavefunction and \ r (R : t) the conditional nu¬ 
clear wavefunction that parametrically depends on the 
electronic configuration r and satisfies the partial nor¬ 
malization condition f dR|;y r (R. t)| 2 = 1 for every r 
at each t. The electronic wavefunction can be shown 
to yield the exact N e — body electronic density, and N e - 
body electronic current-density of the system. The equa¬ 
tions of motion that the electronic and nuclear factors 
satisfy are presented in [ 18 ], The electronic equation, 
in particular, has the appealing form of a TDSE that 
contains an exact time-dependent scalar potential as 
well as a time-dependent vector potential: for any one¬ 
dimensional case, we can choose a gauge such that the 
vector potential is zero 11611181 . and then: 

(-^^2 +e e ( ^,i)) ®(z,t) = id t ®(z,t). (2) 

Hence, the exact potential driving the electron dynam¬ 
ics is e e (z, t ) which can be compared with the traditional 
potentials used to study electronic dynamics. This po¬ 
tential, termed the time-dependent potential energy sur¬ 
face for electrons, e-TDPES, can be decomposed to 

e e (z,t) =e*™(z,t)+T n (z,t) + K.i° nd (z,t) + ef(z,t), ( 3 ) 

where for our H.j" model, 

e app (z,t) = (Xz(R,t) I Wen(z,R) + Wnn(R) \ Xz(R,t)) R 
+ V l (z,t), (4) 

is an approximate potential that generalizes the tradi¬ 
tional quasistatic potential to the case of a quantum nu¬ 
clear wavepacket. The second term, 

T n (z,t) = - (Xz(R,t)\d 2 R \xz(R,t)) R /M , (5) 
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FIG. 2. The exact electronic potential e e (black solid line) and 
its various components together with the quasistatic potential 
e qs (blue dash line) are plotted in the left-hand side at five dif¬ 
ferent snapshots of time. The exact electron density together 
with the electron density calculated from propagating the elec¬ 
tron on e qs and e app are plotted on the right-hand side. 

represents a nuclear-kinetic contribution to the elec¬ 
tronic potential from the conditional nuclear wavefunc- 
tion, while, 

>C c ° nd (z,t ) = (d zXz (R,t)\d zXz (R,t)) R /M , (6) 

represents an electronic-kinetic-like contribution from 
the conditional nuclear wavefunction, and, 

= {xz(R,t)\-id t \xz(R,t)) R , (7) 

is the gauge-dependent component of the potential. It is 
important to note that e app reduces to the quasistatic po¬ 
tential when the nuclear density is approximated classi¬ 
cally as a ^-independent delta-function at R(t) = ( R)(t ): 

e^%z,t\R(t)) = W en (z,R(t)) + W nn (R(t))+V l (z,t). (8) 

We now investigate the e-TDPES (Eqj3| and discuss 
the impact of its components (Eqs. Hfj7T on the elec¬ 
tron dynamics. In particular, we ask how well elec¬ 
tron propagation on e app performs: is accounting for 
the width of the nuclear wavepacket, and its correla¬ 
tion with the electron dynamics via the parametric de¬ 
pendence, enough to capture accurately the full elec¬ 
tron dynamics? In Figure[2]the exact e-TDPES, e e , (black 
solid line) and its various components together with the 


quasistatic potential e qs (blue dash line) are plotted on 
the left-hand side at five different snapshots of time in 
which the field is at the maximum of the cycle. The ex¬ 
act electron density together with the electron density 
calculated from propagating the electron on e qs and e app 
are plotted on the right-hand side. In all these calcula¬ 
tions, the initial state is chosen to be the exact ground 
state density of the complete system. We also plot the 
ionization yields calculated from propagating the elec¬ 
tron on different components of the exact electronic po¬ 
tential on the left on Fig. [3] 

We have chosen times that are representative of three 
different phases of the dynamics (refer to Fig. 0: (1) 
the initial phase, up to t « 20 T, for which the disso¬ 
ciation and ionization probabilities are still negligible 
(panel (a) of Fig. [2j shown at t = 14.25 T), (2) the sec¬ 
ond phase, ~ 20 T < t <~ 35 T, which is when ioniza- 
ton/dissociation mostly occurs; panels (b)-(d) of Fig. [2] 
(t = 24.25, 29.25, 33.25 T) are chosen to illustrate this 
phase, (3) the final phase, t > 35 T, in which the system 
begins to stabilize, while the field intensity decreases, 
represented at t = 36.25 T in panel (e). 

In the first phase of the dynamics (representative 
panel (a)), the nuclear wavepacket is quite localized 
around its initial position and the e-TDPES, quasistatic 
potential and approximate potential are essentially on 
top of each other in the central region (|z| < 10 a.u.). 
They differ from each other only in the tail of the elec¬ 
tronic density, where, in particular, e e has a large step 
downward (Fig. 0a)). Since the density is very tiny in 
the tail region, the overall dynamics is not affected sig¬ 
nificantly by this feature. In the second phase of the 
dynamics the nuclear motion begins to pick up, affect¬ 
ing the shape of the exact e-TDPES in the central re¬ 
gion. From this point on the exact potential begins to 
develop features that are absent in the quasistatic po¬ 
tential. As part of the nuclear density begins to stretch 
apart, the e-TDPES begins to exhibit a double well struc¬ 
ture in the up-field side of the potential (0 < z < 5 
a.u.), while the down-field side maintains a single well 
structure as is shown in panels (b) and (c). Further, the 
depth of the central wells are decreased compared to the 
quasistatic picture. Outside the central region (|z| > 5 
a.u.) the e-TDPES drops down, yielding a barrier that is 
smaller and narrower than that of the quasistatic poten¬ 
tial. This feature in particular, significantly facilitates the 
tunnelling ionization of the electron density in the exact 
dynamics already at t = 24.25 T, evident in the spread¬ 
ing of the exact density (panel (g), see also left panel of 
Fig. |3). In the case of the quasistatic and approximate 
potentials, the ionization is still negligible at this time, 
due to small tunnelling probability. The differences be¬ 
tween the exact e-TDPES and both e app and e qs continue 
to grow in the central region (|z| < 5 a.u.) through¬ 
out the second phase (panels (b)-(d), and correspond¬ 
ing electronic densities (g)-(i)), as contributions from ef d 
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and /Cg° nd get larger and extend closer to the center. It is 
interesting to note that ef d typically has large steps that 
lowers the potential on both sides, allowing for more 
ionization (see also Fig. [ 5 } left panel), while /C'; orid devel¬ 
ops several (smaller) barrier structures, whose net effect 
also appears to increase the ionization yield in this phase 
(see Fig. [3]). The T n term has very small barriers in the 
outer region whose tendency is to confine the density, 
leading to a decrease in the ionization probability. 

By the end of the second phase, at t = 33.25 T (Fig. [2] 
(d)), the exact potential is totally different from the qua¬ 
sistatic potential, everywhere except for at small M < !/ 
presenting a shallow double well structure in both up- 
field and down-field sides of the potential. Furthermore, 
the discrepancy between the e qs and e app becomes more 
noticeable as the nuclear wavepacket begins to split and 
dissociate in the field. By this time, there has been sig¬ 
nificant ionization in all three cases (left panel of Fig. |3}, 
although more in the exact case, as discussed above. To¬ 
wards the end of the second phase, the ionization yields 
of the quasistatic and approximate calculations begin to 
differ from each other, as expected from the growing dis¬ 
crepancy between their respective potentials. 

Entering now the third phase of the dynamics, four 
optical cycles later (panel (e) in Fig. |2}, the exact poten¬ 
tial differs dramatically from the other two showing a 
formation of four wells in the central region (1*1 < 6 
a.u.). The two wells in the center are associated with the 
part of the nuclear density that is not dissociated while 
the other two are associated with the dissociating frag¬ 
ment, and hence, they move outwards. The e-TDPES 
consequently localizes the electronic density in three po¬ 
sitions as seen in panel (j) of Fig. [2| namely in the center 
corresponding to the part of the nuclear density that is 
localized around the equilibrium, and on each of the dis¬ 
sociating fragments of protons. In the third phase, e app 
grossly overionizes the system; as e app has many shal¬ 
low barriers and continues to oscillate in the field, failing 
to stabilize. The quasistatic potential e qs retains a deep 
double well structure throughout the dynamics, in con¬ 
trast to the exact; as the field reduces toward the end of 
the pulse the ionization in either of these cases saturates, 
but the quasistatic fails to get the density and ionization 
yield correct. 

The left panel of Figure [3] shows that neglecting all 
the electron-nuclear correlation terms except for e app un¬ 
derestimates the ionization at first, but later, as the ex¬ 
act ionization begins to saturate, the ionization from 
e app continues to grow, and leads ultimately to a signif¬ 
icant overestimate of the total ionization. Even propa¬ 
gating on the quasistatic potential, a crude approxima¬ 
tion given the earlier discussion, gives a better ioniza¬ 
tion yield. We see from Fig. [3] (left) that adding T n to 
e app reduces the ionization probability at all times, due 
to its small confining barriers as mentioned above. On 
the other hand, adding /C cond to e app increases the ion- 
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FIG. 3. Ionization yields calculated from propagating the elec¬ 
tron on different components of the exact electronic potential 
as well as on the quasistatic potential. Left: A = 800 nm and 
intensity I = 2 x 10 14 W/cm 2 (50-cycle). Right: A = 600 nm 
and intensity / = 10 14 W/cm 2 (20-cycle). Legends apply to 
both. 


ization at first, and then decreases it, giving an over¬ 
all somewhat improved prediction of the ionization dy¬ 
namics relative to dynamics on e app alone. Although 
adding both /C Koml and T n to e app seems to give a good 
final ionization yield, the intermediate dynamics is not 
very good. Adding ef d to e app drastically overshoots 
the ionization, yielding ultimately a completely ionized 
molecule. It appears for the current choice of laser pa¬ 
rameters and initial state all these dynamical electron- 
nuclear correlation terms are important to include to 
obtain good prediction of the ionization yield. But, is 
this conclusion general? Does e app always perform so 
poorly? 

We performed the same calculations using a 20-cycle 
pulse of wavelength A = 600 nm (ui = 0.076 a.u.) and 
intensity / = 10 14 W/ cm 2 , with a sine-squared pulse en¬ 
velope. We further set the initial state to be the 6th ex¬ 
cited vibrational state (c.f. El)- The resulting ionization 
yields computed from the propagation of the electron 
using the different components of the exact potential is 
presented in the right panel of Fig. [3] (T' = 2 fs, is the 
duration of one cycle.). We see that the result of elec¬ 
tron dynamics on e app in this case agrees very well with 
the exact result, and further, that with the addition of 
T n becomes even better. Other combinations of the po¬ 
tential components do not provide satisfactory results. 
In particular, once again ef d severely overestimates the 
ionization, while the the quasistatic dynamics underes¬ 
timates the ionization yield significantly. In this case, the 
exact potential differs substantially from the quasistatic 
potential from the very start, due to the vibrational exci¬ 
tation of the initial state. 

From the discussion above it is clear that ioniza¬ 
tion dynamics depends crucially on the coupling to the 
quantum nuclear motion; accounting for both the split¬ 
ting of the wavepacket as well as its dynamics is im- 
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FIG. 4. Time-resolved, f?-resolved, ionization probability, 
I(R , t). Upper panel: A = 800 nm and intensity I = 2 x 10 14 
W/cm 2 (50-cycle). Lower panel: A = 600 nm and intensity 
I = 10 14 W/cm 2 (20-cycle). 

portant. At this point, one may ask from which part 
of the nuclear wavepacket is the ionization mostly oc¬ 
curring? To answer this, we introduce a time-resolved, 
i?-resolved, ionization probability via 

[ dz\V(z,R,t)\ 2 , (9) 

Jz'i 

with f z , = f ^ taking Z] = 15 a. u., plotted in 

Fig |4] for both of the laser parameters studied in this 
work. In both cases, we observe a clear peak of I(R , t), 
centered around 6 a.u.< II < 7.5 a.u., the region pre¬ 
dicted by the quasistatic analysis of CREI, soon after 
the fields reach their maximum intensities . Hence, 
the quantity I{R, t) represents a very useful measure of 
CREI in a fully dynamical picture, indicating clearly the 
dominant internuclear separations at which ionization 
occurs. We also point out that this quantity is analogous 
to the ionization probability at a given internuclear sep¬ 
aration in the quasistatic picture l24l . 

In summary, we have found the exact potential driv¬ 
ing the electron dynamics in a model one-dimensional 
H 2 molecule undergoing CREI. The presented potential 
provide a complete details of the CREI process beyond 
the quasistatic picture traditionally used to analyze and 
interpret this process. The large differences in the two 
potentials reveals the importance of dynamical electron- 
nuclear correlation terms lacking in previous pictures of 
CREI. The drastic impact these can have on the dynam¬ 
ics has been demonstrated, i.e. propagating the elec¬ 
trons in a potential that neglects these terms lead to large 
errors in the predictions of the ionization yield. Go¬ 
ing beyond the quasistatic treatment by only accounting 
for the width and splitting of the nuclear wavepacket 
is generally not enough to get the correct dynamics of 
CREI. How significant the dynamical electron-nuclear 
effects are for CREI phenomena in larger systems f25l 
remains to be investigated, and moreover, how to ac¬ 
curately model them in practical approximations will 


open a major avenue for future research. We further 
presented a time-resolved, /t-resol ved measure of CREI 
that takes into account the dynamical electron-nuclear 
correlation. This quantity has a clear peak in the re¬ 
gion predicted by the quasistatic analysis of CREI, for 
the fields studied here, soon after the fields reach their 
maximum intensities. 
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